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Abstract 



We study numerically evolution of a self-gravitating nonequilibrium Bose gas 
contained in a fixed volume. We find that, even when the volume is small 
enough to prevent collisionless instability, a compact phase-correlated object 
(a Bose star) can form, by gravity alone, from a disordered initial state. We 
interpret this effect and the associated growth of the density contrast as conse- 
quences of a new type of gravitational instability — a nonlinear instability due 
to stimulated gravitational scattering of the bosons. Our results imply that 
formation of Bose stars, in regions that have fallen out of the Hubble expan- 
sion, may be a quite general phenomenon, not requiring a large preexisting 
correlation length or any short-range interactions. 
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I. INTRODUCTION 



Gravitational (Jeans [|I|) instability is believed to be at the root of all the large-scale 
structure that we observe in the universe. (For a review on structure formation, see ref. 
The instability is a direct consequence of gravity being an attractive force: particles will 
accrue on small clumps, and the clumps will get bigger. Jeans's original calculation was based 
on description of matter via classical hydrodynamics. This description applies only at spatial 
scales larger than the mean-free path of the particles (which is defined here with respect to 
some interaction other than gravity). Evolution of density perturbations with scales smaller 
than the mean-free path can for many purposes be described by the collisionless Vlasov 
approximation — unless the scale is so small that this semiclassical method breaks down. In 
the collisionless approximation, a self-gravitating gas is unstable with respect to sufficiently 
long-wave perturbations, and the instability is of the same type as Jeans's. (We will define 
later more precisely what "the same type" means.) 

For a wide class of initial distributions, however, the collisionless instability disappears for 
short-scale perturbations. Consider, for example, a homogeneous (on average) nonrelativistic 
gas of particles of mass m in a box of fixed size, in the Newtonian approximation to gravity 
(in addition, we assume that the uniform component of density does not gravitate). We show 
in Appendix that, for the class of momentum distributions defined there, the collisionless 
instability exists only for density perturbations with wavenumbers 

k < ckyko = k, , (1) 

where fco is a typical momentum of the particles, and C is a numerical constant depending 
on their momentum distribution. (We use a system of units with ^ = 1.) The wavenumber 
kj is determined by the average density number nave: 

k'^j = 167rGm^nave , (2) 

and will be called the Jeans wavenumber. One may wonder if any gravitational instability 
exists for perturbations with /c > /c*, that is at shorter scales. 

The condition (|I|) suggests that, depending on the relation between ko and kj, a self- 
gravitating gas can be in one of two different regimes. When 

fcj > fco , (3) 

eq. (|1|) shows that collisionless instability occurs for perturbations with all wavenumbers 
up to ~ kf), at which point the semiclassical description of particles breaks down, and 
the Vlasov equation can no longer be used. In at least one case, however, one can show 
that a linear instability exists even for k > ko (although of course the argument cannot 
be based on the Vlasov equation). This case is a perfectly ordered Bose gas, i.e. a gas in 
which all particles are Bose-Einstein condensed in the k = mode. Elementary excitations 
of this Bose-Einstein condensate are Bogoliubov's quasiparticles, and their dispersion law 
is known, from Bogoliubov's work 0, for arbitrary two-body potential. When the only 
interaction between particles is Newtonian gravity, the dispersion law (for a nonrelativistic 
gas) reads 
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uj{k) = (l/2m) [fc^- , (4) 

which means that the perfectly ordered Bose gas is gravitationally unstable for all k < kj. 
This result had been reported previously in a number of papers (We have presented it 

here for a region that is out of the Hubble expansion; a treatment in the expanding universe 
is also available [^H-) 

The perfectly ordered state is just one possible initial state of a Bose gas. A generic initial 
state is a patchwork of correlated domains, rather than one such domain. If we describe 
a nonrelativistic Bose gas by a complex field ip, correlated regions are those in which the 
phase of ip is more or less constant in space. The typical size of correlated regions is 27r//co. 
If such a region contains many particles, we can apply formula @) to it individually, but 
only for wavenumbers k ^ k^: only these sample a more or less ordered background. Eq. 
(1) then predicts an instability for k in the range k^ <^ k < kj. At A; ~ fco; we expect 
this instability to cross over smoothly to the instability found in the Vlasov approximation. 
Thus, we expect that a nonrelativistic Bose gas satisfying the condition (^) is unstable with 
respect to density perturbations with all wavenumbers k < kj (although the theoretical 
description of the instability has to be switched at ~ k^). 

We can now formulate more precisely what we mean by saying that the gravitational 
instability in the perfectly ordered Bose gas is of the same type as the instability found 
by Jeans for matter described by classical hydrodynamics. The similarity has two closely 
related aspects. First, in either case, the instability is linear, i.e. density perturbations with 
different wavenumbers do not have to interact with one another for it to occur. Second, the 
growth rate of the instability, for low-fc modes, scales as square root of the average number 
density: 

lmu;(0) oc ^/n^ (5) 

(see e.g. eq. (§)). One can say that this square- root scaling law defines a certain "universality 
class" , to which the original Jeans instability, the coUisionless instability, and the instability 
of the ordered Bose gas all belong. 

In this paper, we want to consider a nonrelativistic Bose gas in which the relation between 
ko and kj is reversed compared to (H): 

kj<^ko. (6) 

Our motivation is simply that we do not know at present which type of initial state is more 
likely to be realized in the early universe. We will call the case (|^) disordered (because 
correlated domains are initially small), as opposed to the case (|^), which we will call ordered 
(initial correlated domains are large). According to eq. (|^), the coUisionless instability 
will occur in a disordered gas only at large spatial scales. So, if there is a gravitational 
instability on shorter scales, it is caused by gravitational scattering of the particles, rather 
than by coUisionless effects. Our goal was to determine effects of gravitational scattering on 
the evolution of a disordered Bose gas. 

We imagine the following sequence of events in the early universe. A disordered Bose gas 
first becomes unstable with respect to large-scale density perturbations via the coUisionless 
(Vlasov) mechanism. As clumps formed by this mechanism grow denser, coUisionless insta- 
bility can develop on shorter scales. We expect that this process will give rise to clusters of 
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matter, each cluster including a number of smaller clumps and some inter-clump gas. Grav- 
itational scattering, which is a slower process, will operate inside those clusters, after they 
have already fallen out of the Hubble expansion. Accordingly, we considered a disordered 
Bose gas in a box of fixed size. In addition, we used Newtonian gravity, which we expect to 
be good for sufficiently short-scale perturbations. 

Our results, obtained by numerically following the evolution of a random classical field, 
apply when the initial state has large (compared to unity) occupation numbers at low mo- 
menta and little occupation at high momenta. These results are as follows, (i) Starting 
from a disordered state, we have observed a well-defined crossover from coUisionless insta- 
bility to a slower growth of the density contrast. This slower growth can be viewed as a 
gravitational instability of a different type than the linear instabilities discussed above. We 
confirm the existence of this new type of gravitational instability by running simulations 
in smaller boxes, so as to remove the coUisionless instability altogether. In this case, any 
growth of the density contrast will have to be attributed to some new type of instability. We 
indeed observe such a growth, and we interpret it as a result of gravitational scattering of 
the bosons. The scaling of the rate of the growth with riave, in a small box, is quite distinct 
from the square-root law (^). In this sense, the new instability is in a different "universality 
class" than the linear instabilities, (ii) Clumps that result from this nonlinear instability 
are phase-coherent; in other words, the instability leads to formation of Bose stars. (By a 
Bose star we mean any coherent gravitationally bound configuration of a Bose field; ground 
states of gravitating bosons were originally discussed in ref. [0].) It has been long considered 
possible that Bose stars will form as result of the linear instability in an ordered Bose gas 
(although, to our knowledge, that had not been actually demonstrated). Formation of 
Bose stars in an initially disordered gas is a different effect. We attribute it to gravitational 
scattering in our case being stimulated: as an effect of Bose statistics, particles that scatter 
into regions of high density will prefer to end up in phase with particles that are already 
there. 

The rest of the paper is organized as follows. Sect. 2 contains formulation of the problem 
in terms of a random classical field. In Sect. 3 we present numerical results. Our conclusions, 
and connections to earlier work, are given in Sect. 4. Some technical details concerning the 
coUisionless instability are given in Appendix. 



II. SETTING UP A RANDOM CLASSICAL FIELD 

We consider a nonrelativistic Bose gas described by the following equations: 

dih 

i-^ = -—^lj + m^^ + g^i/j^ijj'^ ; (7) 
= 47rGm(V'V - ^ave) • (8) 

Here ip is a. complex Bose field, and $ is the (real) gravitational potential. Number density 
of bosons is n = ip'^ip, and its volume average is called riave- For reasons given in the 
introduction, we assume it sufficient to consider the gas in a box of fixed size (rather than 
in an expanding universe) and to use the Newtonian approximation to gravity. 

All numerical results presented in the next section are for = 0, i.e. for the case when 
gravity is the only interaction between the particles. The method that we use can be applied 
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equally well to the case ^ 0. In this paper, however, we limit ourselves to only a cursory 
discussion (in the concluding section) of possible effects of short-range interactions. 

The field ip is in principle a quantum operator but here we treat it as a random classical 
field. Usefulness of this classical approximation has been recently emphasized in connection 
with a number of problems in nonequilibrium statistics of Bose fields. Among these problems, 
the one closest to the problem at hand is phase separation in a nonequilibrium Bose gas 
with a local attractive interaction (and a repulsive interaction at larger densities) |^ . Some 
of the results that we obtain here, notably, the coherent nature of the clumps, are similar 
to those for a local interaction, but due to the long-range nature of gravity there are also 
some differences. For example, for a short-range interaction, one can also identify ordered 
and disordered types of initial states. Results of ref. correspond to a disordered initial 
state. For neither type of state, however, we expect a large-distance crossover to collisionless 
(Vlasov) dynamics. 

The system (0)-(|^) can be viewed as a modification, required to include Newtonian 
gravity, of the Gross-Pitaevskii (GP) equation familiar from the theory of laboratory Bose 
gases |T^. The system (0)-(H) is universal in the same sense as the GP equation is: at 
low gas densities, the only parameter needed to characterize short-range interactions is the 
scattering length, which is encoded in the coupling constant g^. 

The classical approximation is good when the field is large and its power spectrum 
is contained at small momenta: speaking in terms of particles, the occupation numbers in 
low-momentum modes are much larger than unity, while those in high-momentum modes 
are small and rapidly decrease with momentum. We stress, however, that when scattering 
is strong, as it may become at later stages of clumping, the very notion of a particle, as an 
entity propagating for relatively large times without collisions, becomes ill-defined. One of 
the useful features of the classical approximation is that it allows one to avoid dealing with 
that notion altogether. In addition, provided the power spectrum of ip satisfies the necessary 
conditions, the classical approximation will apply equally well either in the ordered or in 
the disordered case. We also stress that, although this may sound paradoxical, the classical 
approximation does take into account effects of the quantum Bose statistics of the particles. 
Indeed, the very possibility for ip to become (almost) classical, in the limit of large occupation 
numbers, is one such effect. 

For numerical simulations presented in this paper, we have chosen the initial power 
spectrum of ip in the following simple form (cf. ref. 0): 

|^kl' = ^exp(-A;VA:o'), (9) 
where ipy^ are Fourier components of the field ip: 

^(r) = \/-V2^^j^exp(zkr), (10) 
k 

V is the integration volume. For an ideal gas, iV^j^P would be the conventionally defined 
occupation numbers; we will often use the same terminology for the interacting gas. For the 
classical approximation to apply, it is sufficient that A ^ 1. Although the distribution (^ 
may look like an equilibrium Maxwell distribution, it is totally unrelated to the latter. The 
Maxwell distribution follows from the Bose distribution in the limit when the occupation 
numbers are small. In our case, A ^ 1, and the occupation numbers (at low wavenumbers) 
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are large. Thus, the spectrum (j^) is highly nonthermal. Eq. determines the magnitudes 
of The phases of are chosen as uncorrelated random numbers. The initial number 
density contrast corresponding to these initial conditions is 5n/n ~ 1. 

Bose gases with highly nonthermal spectra can form in the early universe through vari- 
ous nonequilibrium processes, such as spinodal decomposition or postinflationary reheating. 
(Formation of nonthermal spectra in the latter case is reviewed in We do not pretend 

that the specific form (^) is in any way realistic. We only use it as a representative of 
spectra satisfying the conditions of applicability of the classical approximation. We have 
also considered some different forms of the initial spectra and have convinced ourselves that 
these modifications do not alter the results in any significant way. 

Finally, for comparison with the classical method used here, we review results obtained in 
the collisionless (Vlasov) approximation. Evolution of density perturbations with wavenum- 
bers k <^kQ can be studied by considering particles as semiclassical wave packets. We stress 
that this semiclassical approximation for the particles is entirely different from the classical 
approximation for the field if). In particular, it does not require that the occupation numbers 
be large but requires that perturbations be slowly varying. (It is the other way around for 
the approach based on a random classical field.) When the semiclassical view of the parti- 
cles applies (i.e. a perturbation is slowly varying), one can describe them with a classical 
distribution function /(r, p; t), which is proportional to the Fourier transform, with respect 
to k, of the quantum average of ^jj'^ k/2^P+k/2' "^^^ distribution function, together with 
the gravitational potential $, can be used to construct a collisionless approximation (which, 
if needed, can be augmented by a collision integral). The construction is entirely analogous 
to the derivation of the Vlasov equation for plasmas, see e.g. |jl2|. In the collisionless ap- 



proximation, we find that frequencies ci;(k) of small perturbations of a homogeneous state 
with a distribution function /o(p) are determined by the dispersion equation 

i + l^/k§^^!!^ = o. (11) 

J ap kv — — zO 

where v = p/m. When kv, for a typical f , is much smaller than \uj\, the dispersion equation 
reduces to 

1 + 47rG'mnave/c^^ = . (12) 

Eq. ([12D shows that there is a linear instability with the growth rate 

cJi(O) = (47rG'mnave)^/^ = k]/2m . (13) 

For a wide class of initial distribution functions, which is specified in Appendix, the colli- 
sionless instability disappears at /c > fc*, where fc* is estimated as in the condition (|l|). The 
precise value for k^, is obtained by setting = in eq. ( pA] ) and solving for k. (This proce- 
dure is justified in the Appendix.) For example, in an initial state with power spectrum (^ 
the distribution function is 

/o(p) = exp(-pVPo) , (14) 

where po = in the system of units with h = 1. In this case, collisionless instability occurs 
only for 
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k<K = {^/2y/^ko , (15) 
where ^ = kj/k^ Eq. (|15D is equivalent to (|) with C = 1/^2. 

III. NUMERICAL RESULTS 

For numerical work, it is convenient to rescale the fields of our model, so that we have 
to deal with a smaller number of parameters. Define new fields x ^i-nd as 

X = (87rGm3)i/2^ , (16) 
(f) = 2m^^ . (17) 

Then, eqs. (0)-(§) with = take the form 

2rm^ = -V\ + <Px , (18) 

= X^X - '^ave , (19) 

where z/ave is the average oi u = x}x over space. Initial condition (^ takes the form 

|XkP = Sexp(-A;V^o) , (20) 

where B = SnGm^A, and Xk are related to in the same way as ip]^ to ip{r), see eq. 
(0). Finally, the expression (0) for the Jeans wavenumber reduces to 

k^J = 2j/ave . (21) 

All data presented in the figures were obtained by numerically integrating the system 
(p!6|)-(p^), with initial power spectra of the form (^) and uncorrelated random initial phases 
for on 64'^ cubic lattices with periodic boundary conditions. We have chosen the unit 
of length so that fco = 2n, and the unit of time so that 2m = 1. The integrations were 
done using a second-order in time operator-splitting algorithm; updates corresponding to 
the operator were carried out by a spectral method based on the fast Fourier transform. 
The algorithm conserves the number of particles exactly; energy nonconservation was below 
1% for all the data sets represented in the figures. 

As we have mentioned in the introduction, there are basically two choices for the size of 
the integration box. One possibility is to choose the box large enough to activate collisionless 
instability, in order to see a crossover from the collisionless instability to a different regime. 
While it is clear, even a priori, that a linear instability cannot carry on forever and will 
have to end somehow, it is still interesting to see explicitly how that happens and what the 
different regime is. Representative results are shown in Fig. 1. These results are for B = 20 
and box size L = 4.25. A crossover, around t = 1.2, from a more rapid to a slower overall 
growth is seen both in the density contrast and, especially, in the ratio of the total potential 
and kinetic energies. The density contrast is defined as 

6n/n={{u-u,,,)YV^..o, (22) 
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where the brackets denote averaging over space. It is significant that the growth of the 
density contrast continues after the hnear instabihty apparently terminates. This allows 
us to speak about a nonlinear instability, which we attribute to gravitational scattering. 
The overall growth of the density contrast is superimposed on rapid oscillations, which are 
especially large at the later, scattering stage. We attribute these oscillations in the density 
contrast to oscillations of matter about a forming star. 

To confirm the existence of a nonlinear instability, it is useful to run simulations in 
a smaller box, so that the long-wave modes that could undergo the collisionless (linear) 
instability are eliminated. For a cubic box with periodic boundary conditions, this means 
choosing the side length L so that 27r/L is larger than the wavenumber k^, that appears 
on the right-hand side of (|l|) (there are no density perturbations with k = because the 
total number of particles is conserved). Any growth of the density contrast in an initially 
disordered gas in such a box will have to be attributed to some new type of gravitational 
instability. We can find out if this instability belongs to a new "universality class" by 
studying how its rate depends on the initial density of the gas. 

When gravity is the only interaction between the particles, the rate of clumping, during 
an initial stage when the gas is still nearly uniform, can be written in general as 

t-' = §^F{kykl k.L) . (23) 

The function F depends on two parameters that measure the strength of gravitational 
interaction between the particles: 

e = k)lkt (24) 

measures the strength of interactions for momentum transfers of order 27r/fco, while k^L 
measures that for momentum transfers of order k^am = 2tt/L. In the limit k^L oo, 
when the initial clumping is due to collisionless instability, F becomes a function of its first 
argument only, and 



no « ve , (25) 



see (|T2D. In an ordered Bose gas, dependence of the rate on L disappears whenever kj 3> 
/cmin, and it scaling with ^ is then given by the same square-root law (|23|), see (^). 

In an initially disordered gas, situation is more complicated. The condition (^ (i.e. 
^ 1) guarantees only that interactions in the initial state are weak for momentum transfers 
of order 27r//co. In a small box, we also have the condition k^, < /cmiru which guarantees that 
interactions are at most of moderate strength even for the smallest momentum transfers. 
In this case, the initial gas can be considered weakly interacting, and one may contemplate 
applying kinetic theory. However, even if suitable kinetic equations exist, they will, in 
general, not reduce to a Boltzmann equation (properly modified to include the effect of 
large occupation numbers). Indeed, as we have already noted, fluctuations with spatial scales 
£ ^ 27r/A;o cannot be described using the conventional semiclassical distribution function. If 
such fluctuations are important (and we will argue that in the present case they are), the 
Boltzmann equation will not apply. 

The simplest way to see that the Boltzmann equation is in general insufficient for studying 
evolution of a Bose gas, either for long-range or for short-range interactions, is to note that it 
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includes interaction only via the scattering probability and so does not distinguish between 
attractive and repulsive interactions. On the other hand, on physical grounds we expect 
dynamics for these two types of interactions to be vastly different. 

Let us refer to the rate predicted by the Boltzmann equation for changes in the dis- 
tribution function as the Boltzmann rate. For a short-range attractive interaction, it has 
been found numerically in ref. that scaling of the clumping rate with ^ agrees to a good 
accuracy with scaling of the Boltzmann rate. In general, however, there is no reason to 
expect such agreement. For gravitational scattering, the Boltzmann rate (for a gas with 
large occupation numbers) is of the form (PB|) with 



F{^,KL) (x^Hn{koL) . (26) 

The logarithm here is analogous to the Coulomb logarithm in plasma kinetics (cf. ref. [0); 
it is due to an infrared divergence in the collision integral. In a small box, this divergence 
is cut off at wavenumber transfers of order fcmin = 2ti/L. Now, consider the case when 
and L are fixed. The scaling predicted by ( PB]) is F{^) oc This is not quite what we see 
numerically. Instead, a good fit to the rate is obtained with 

F(0 (X r • (27) 

where a is noticeably larger than 2. 

We attribute this large deviation of a from 2 to an important role played by fluctuations 
with spatial scales much smaller than l-njkQ. One can readily invent a reason why such 
fluctuations are more important in the presence of a long-range interactions than in the 
absence of such: even a short-scale fluctuation will influence many particles when there is a 
long-range force. We thus arrive at a picture of short-scale fluctuations constantly appearing 
and disappearing, until one of them is able to start a growing clump. 

To extract the scaling law from our numerical results, we have used data corresponding 
to L = 3.5 and three different values of 5: B = 7, B = 10 and B = 20. The corresponding 
values of ^ = fcj/^o 0.050, 0.071, and 0.143. For these values of L and ^ , the right-hand 
side of (|15D is smaller than 27r/L, so coUisionless instability does not occur. Consequently, 
the density contrast's growth, which was observed in all three cases, is a signature of a new 
type of gravitational instability. 

The growth of the density contrast should be more accurately referred to as a growth of 
some average with respect to time, because it is superimposed on rapid oscillations similar 
to those seen at later stages in Fig. 1. We expect that such an average will depend on 
time and ^ mainly through a dependence on the ratio t/tc, or, because k^ and L are kept 
constant, through the product t^", where a is the power that we want to determine: 



6n/n ^ f{tC) . (28) 

The overline denotes averaging over several oscillations, and / is some function. When we 
plot the averaged density contrast, for the above three values of B, against the rescaled time 
variable t^", we find that a good collapse of the plots is achieved for 

a = 2.7 , (29) 
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see Fig. 2. The collapse is not nearly as good for a = 2, the value characteristic of 
a Boltzmann rate. Our interpretation is that the instability is a result of gravitational 
scattering, in which fluctuations with spatial scales £ -C 27r//co play an important role. 

There is an approximately linear stage in the growth of 5n/n in Fig. 2, from about 0.5 
to about 1.2 in units of lOOt.^^'^. Fitting this linear growth with a time dependence of the 
form 



5n/n = t/tc + const , (30) 

provides one possible definition of the rate t~^. Numerical fitting gives 

t-^ ^ 250^2-^ . (31) 

This corresponds to eq. (p^) with -F(^) ^ 6^^ ''. 

The estimate (|3lD has been obtained for initial power spectra of the form (P). It is 
natural to expect that a different form of the initial power spectrum will lead to a different 



numerical coefficient in (|3T|). That is indeed so, although in cases that we have considered 
the difference is not overwhelming. For instance, one can generate a random field with a 
power spectrum of the form (P) and then make, by hand, the magnitudes of x on all lattice 
sites equal to some fixed value, while preserving the phases (in coordinate space). One can 
then use this new field as an initial condition that has a power spectrum different from (^. 
We have done that for L = 3.5, and i? = 10 and 20. In this case, the density contrast is 
initially zero. However, it rapidly develops to a nonzero value, as dynamics of the density is 
activated by dynamics of the phase of ip. The subsequent evolution of the density contrast 
is fairly similar to that seen in Fig. 2. The scaling of the rate is consistent with a = 2.7, 
but the rate itself is about half of that given by eq. (|3T|). 

Finally, Fig. 3 shows the profile of the star at different moments of time, and Fig. 4 
illustrates the star's phase-correlated (Bose-Einstein condensed) nature. (Only one large 
star has formed per box, in addition to a number of much smaller clumps.) 



IV. DISCUSSION 



The main result of this work is numerical evidence for a new type of gravitational insta- 
bility. We have observed this instability in simulations of a disordered Bose gas in a box of a 
fixed size. We interpret it as being due to collisions between the particles, or, in terminology 
more suitable for a system with large occupation numbers, between classical waves. We have 
found that scaling of the instability rate with the average density deviates from the behavior 
characteristic of changes described the Boltzmann equation. We interpret this deviation as 
an indicator of an important role played by short-scale fluctuations of the field. We have also 
presented evidence (Fig. 1) that this instability will take over from a coUisionless instability 
in a large system. This means that it could contribute to formation of structures inside 
regions that had already fallen out of the Hubble expansion. In addition, we have found 
that this coUisional instability leads to formation of a phase-correlated clump of matter (a 
Bose star) in an initially disordered Bose gas. Hence, Bose stars can form through gravity 
alone even in the absence of a large preexisting correlation length. 
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These results were obtained using the classical approximation for a Bose field. The clas- 
sical approximation applies to a field that has a certain type of nonthermal power spectrum, 
namely, occupation numbers are large at small wavenumbers but rapidly decrease towards 
the ultraviolet. We do not exclude, however, that a similar nonlinear instability exists in 
matter with different distributions of particles with respect to momenta, for example, in a 
thermal or nearly thermal gas (for which the classical approximation does not apply). 

Our results may apply to nonbaryonic dark matter (if the dark-matter particle is a boson) 
and to low-density hydrogen. (At high densities, when inelastic processes are important, the 
simple description ([7p is no longer adequate.) It is fascinating to think that some of the 
matter in the universe may be in the form of superfiuid hydrogen. However, there may be 
interesting implications for structure formation even if our results apply only during initial 
stages of clumping, when the density is still low. 

We should also comment on a possible role of short-range interactions in a disordered 
self-gravitating Bose gas. As we have already mentioned, at low enough densities all the 
short-range interactions will be encoded in a single coefficient (74 of the cubic term in eq. 
(0). The strength of this local interaction relative to the strength of gravitational scattering 
in the initial state is measured by the parameter rj = gJ^Q/ ATrGim? . This parameter can be 
small even when (74 is much larger than G/c^ (c is the speed of light), provided the ratio 
ko/mc is small enough, i.e. the initial state is sufficiently nonrelativistic. On the other hand, 
a not-too-small 1] can lead to interesting effects. A repulsive interaction ((74 > 0) in a Bose 



gas causes coarsening — growth of the sizes of correlated domains [|T3|. If the interaction is 



strong enough, we expect that the gas will become sufficiently ordered before the nonlinear 
instability had a chance to develop, and clumping will take place via the corresponding 
linear instability. In addition, we expect that repulsion will make emerging Bose stars larger 
in size and smaller in density. (For equilibrium configurations, it has been long known that 
even a small repulsion makes a large difference |l^.) An attractive interaction {g^ < 0) 
produces correlated clumps all by itself 0; we expect that it will help gravity along when 
both are present. 

Finally, we discuss relation of our work to earlier work in the literature. We can discern 
two trends in the earlier work on formation of Bose stars. Some researchers considered 
gravitational stability and evolution of already coherent configurations. That work included 
linear stability analysis of a homogeneous Bose-Einstein condensate and a study of 
spherically symmetric collapse of an already coherent spherically symmetric configuration 
|T^. It has been suggested [^,|| that the Jeans- type instability of a homogeneous Bose- 
Einstein condensate may lead to formation of Bose stars. Although we do not present 
corresponding data in this paper, we have confirmed that a Bose star indeed forms in an 
ordered Bose gas. However, we have given evidence that a large preexisting correlation 
length, i.e. a high initial degree of spatial coherence, is not necessary for Bose star formation. 
A Bose star can form even in an initially incoherent (disordered) gas via a new type of 
gravitational instability that we have identified here. The second trend in the earlier work 
was to consider ordering effects of local interactions. In particular, it has been argued that 
due to effects of the Bose statistics even a very weak interaction, such as that between 



axions, can result in a relaxation time comparable to the age of the universe |T6|, and, as 



a consequence, a Bose star may form out of an initially incoherent clump |T^. Here we 



have shown that phase coherence can develop due to gravity alone, in the absence of any 
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additional interaction. 
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APPENDIX: MORE ON THE COLLISIONLESS INSTABILITY 

In this appendix we will show that for a wide class of initial distribution functions the 
collisionless instability (such as found in the Vlasov approximation) disappears at sufficiently 
large wavenumbers. Whether the instability is present or absent is determined from the 
dispersion relation (0), which we rewrite here as 

7(^,fc) = 0, (Al) 

where 7 is the "gravitational permittivity" : 

, , , ^ iirGm^ , q'iPx) / * „x 

^{u;,k) = l + —-^ dp^ "^^y . A2 

/c^ j-00 Px — muj/k — t(j 

We have oriented coordinate axes so that k = (/c, 0, 0), with k > 0, and defined a distribution 
function with respect to p^- 

g{Px) = J dpydpJo{p^,Py,p,) . (A3) 

In an isotropic medium, the form of g{px) would not depend on the direction of x (that is 
the direction of k), but we do not assume isotropy here. 

By its physical meaning, g{px) is nonnegative. We assume that for any direction of k the 
function g{px) satisfies the following conditions: (i) g{px) is smooth (infinitely different iable) 
and decreases with \px\ rapidly enough for the integral in ( |A^ ) to be convergent for any lo 
with Imuj > 0. (ii) g{px) is even; (iii) g{px) is monotonically decreasing for all Px > 0. 
From these conditions, two more follow: (iv) g'{0) = 0; and (v) g"{0) < 0. Primes denote 
derivatives with respect to p^- 

Eq. ( [A^ ) defines an analytic function in the upper half-plane of complex u; in our 
definition, the upper half-plane does not include the real axis. To obtain 7(0;, k) in the lower 
half-plane, we have to analytically continue it there from the upper half-plane. (Using eq. 
( |A^ ) directly in the lower half-plane would give another, unphysical sheet of '-f{uj, k).) 

Because g{px) is even, g'{px) is odd. From this, it follows that in the upper half-plane 
^{uj,k) can have zeroes only on the imaginary axis. For u = iuji with cjj > 0, 7 takes the 
form 

,(..a)^l + 4.G™V~*.p|f|L (A4) 

(in particular, it is purely real). The conditions on g{px) guarantee that Pxg'iPx) < for 
any nonzero Px- That means that the integral in ([A4|) is negative, so that at any fixed A; > 
'y{iuJi, k) is a monotonically increasing function of tUj. Consequently, at a given k it has at 
most one zero in the upper half-plane of a;. This zero, when it occurs, signals an instability 
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of the initial state with respect to fluctuations with that k. For example, such a zero always 
exists for /c — i> 0; its position, ^^(O), is given by eq. (p^Sl) . 

Now let us see what happens as k increases. Zeroes of 'y{u,k), i.e. roots of eq. ([Al[), 
define a dispersion law uj{k). As long as, for a given k, the zero is still in the upper half-plane, 
we can continue to use eq. ([A^). First, it follows from (Kl) that ci;i(0) is an upper bound on 



iUi{k). So, as we change k, a root cannot appear from infinity. Second, differentiating (|A1|) 



with respect to k and using again (|A4D, we find that, when g{px) obeys the stated conditions. 



duji{k)/dk is necessarily negative. So, as k increases, the root moves down towards the real 
axis. It reaches the real axis at the value k = k^ determined by setting cuj = in (|A4]): 



ATiGrrf 



(A5) 



(Note that there is no divergence at = here, because g\Q) = 0, see condition (iv).) We 
have already established that any root of (|A1|) in the upper half-plane is on the imaginary 
axis, and that such a root cannot appear from infinity. As a consequence, the point u; = is 
the only point where the root of ([Al|) can leave or enter the upper half-plane of u. Because 



the corresponding value of k is unique, the root will not reappear in the upper half-plane at 
any k > k^,. We conclude that the coUisionless instability is absent for all k > k^,. 

A detailed picture of how the instability disappears can be obtained by expanding the 
original expression for 7, eq. (|A^), in an asymptotic series about uj = 0. The iO prescription 
in (|A2|) will ensure that we are expanding on the correct sheet of 7(ci;, A;): 

AnGm^ , q'iPx) ( muj/k ^, n.\ . 



k'^ J-Qo Px — iO \ Px — iO 
For a smooth even g{px), this evaluates to 

Writing k = k^ + Ak and expanding in small Ak, we find that the root of the dispersion 
equation has the form 



u;{Ak) = -I- 



IZodpxg'{Px)/Px 



9"{0) 



0{{Akf) . (A8) 



This expression shows explicitly that, as k increases past fc*, the root moves from the upper 
half-plane into the lower half-plane. 

The structure of higher-order terms in ([A7|) is as follows: the coefficients of even powers 
of u are all real, while the coefficients of odd powers are all imaginary. So, to any finite 
order in Ak, the root of the dispersion relation will remain purely imaginary even after it 
moves in the lower half-plane, and the mode corresponding to that Ak will be overdamped 
(non-oscillatory). This conclusion may be obviated by terms that are "nonperturbative" in 
Ak, i.e. are not seen in the asymptotic expansion {\KBi). We also stress that it apphes only 
to fluctuations about a homogeneous state, not to a state that already contains clumps. 
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FIGURES 
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FIG. 1. Illustrates crossover from the collisionlcss instability to a scattering regime. These 
plots are for B = 20 and L = 4.25. Solid line shows the density contrast; dashed line — the ratio of 
the total potential and kinetic energies. The crossover occurs around t = 1.2. 
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to 




FIG. 2. Density contrast averaged over intervals of time, as a function of rescaled time. Each 
symbol is at the center of an interval over which the averaging is carried out. All data are for 
L = 3.5. The data for B = 7 and B = 20 correspond to the same realization of random initial 
phases; those for B = 10 — to a different one. 
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FIG. 3. Profile of \iJj\ along a line passing through the maximum of {ipl, at different moment 
of time. Because the location of the maximum moves, the actual lines are different at different 
moments (although they are chosen parallel). These plots are for L = 3.5 and B = 7. The 
coordinate along the line is in units of the lattice spacing. 
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FIG. 4. Distribution of the field ■0 on a quarter of a spatial slice at time t = 100, for L = 3.5 
and B = 7. The slice cuts through the star's center. The length of an arrow corresponds to |0| on 
a given site, and the angle the arrow makes clockwise from 12 noon — to the phase of ip. Arrows 
on all sites with \il>\ > 70 have the same (maximal) length. 
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